function [out] = efe(para, reg)

% Returns [3x2] matrix
%
% e.g. for reg = 1
% out = [efe(y1MSNK_a-y1F_a) efe(pi1MSNK_a-pi1F_a);    /Technology shock
%        efe(y1MSNK_u-y1F_u) efe(pi1MSNK_u-pi1F_u);    /Markup shock
%        efe(y1MSNK_e-y1F_e) efe(pi1MSNK_e-pi1F_e)]    /MP shock

alpha1 = para(1);
alpha2 = para(2);
gamma1 = para(3);
gamma2 = para(4);
beta1  = para(5);
kappa1 = para(6);
tau1   = para(7);
rho_a  = para(8);
rho_f  = para(9);
rho_i  = para(10);
p11 = para(11);
p22 = para(12);


%----------------------------------------------------------------------
% MSNK solution
%----------------------------------------------------------------------
A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
            -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];           
        
solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
            -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];         
solvec_f = inv(A_Af_mat) * [0;0;-kappa1;-kappa1;0;0];

A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
            -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];        
solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];

y1a = solvec_a(1);
y2a = solvec_a(2);
pi1a = solvec_a(3);
pi2a = solvec_a(4);
y1f = solvec_f(1);
y2f = solvec_f(2);
pi1f = solvec_f(3);
pi2f = solvec_f(4);
y1i = solvec_i(1);
y2i = solvec_i(2);
pi1i = solvec_i(3);
pi2i = solvec_i(4);

%----------------------------------------------------------------------
% Fixed-regime solution
%----------------------------------------------------------------------
p11 = 1;
p22 = 1;

A_Aa_mat2 = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
            -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];           
        
solvec_a2 = inv(A_Aa_mat2) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

A_Af_mat2 = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
            -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];         
solvec_f2 = inv(A_Af_mat2) * [0;0;-kappa1;-kappa1;0;0];

A_Ai_mat2 = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
            -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];        
solvec_i2 = inv(A_Ai_mat2) * [0;0;0;0;1;1];

y1a2 = solvec_a2(1);
y2a2 = solvec_a2(2);
pi1a2 = solvec_a2(3);
pi2a2 = solvec_a2(4);
y1f2 = solvec_f2(1);
y2f2 = solvec_f2(2);
pi1f2 = solvec_f2(3);
pi2f2 = solvec_f2(4);
y1i2 = solvec_i2(1);
y2i2 = solvec_i2(2);
pi1i2 = solvec_i2(3);
pi2i2 = solvec_i2(4);

if reg==1   
    out = [y1a-y1a2 pi1a-pi1a2; y1f-y1f2 pi1f-pi1f2; y1i-y1i2 pi1i-pi1i2];
elseif reg==2
    out = [y2a-y2a2 pi2a-pi2a2; y2f-y2f2 pi2f-pi2f2; y2i-y2i2 pi2i-pi2i2];    
end





